
% MATLAB script
% Plot of photodisruption diagram - uses field of simulations to interpolate 
% in logspace onto specific thresholds (boiling, melting, etc.)
% SWH 2011

close all
clear all

% specific case parameters
ro = 25e-9;
Qabs = 0.15;
taup = (200e-15)*(2/log(2))^(1/2);
Fnp = 2534*(0.2274/0.4677)^2;  % adding scaling from approximate corrected results using distribution
Fnpe = 274;
bc = Qabs*Fnp*pi*(25e-9)^2/(70*300^2*4/3*pi*(25e-9)^3);
gc = 2e16*(200e-15)*taup/(70*300);

% constants
g = 2.00E+16;
G = 1.05E+08;
To = 300;
Ceo = 70*To;

% read in of T solution & parameters
sim_min = 1001;
sim_max = 1500;
sim_tot = 500;  % formerly yytot
radii = load('..//inputs//radii.txt');
fluences = load('..//inputs//fluences.txt');
pulses = load('..//inputs//pulses.txt');
rlength = length(radii);
flength = length(fluences);
plength = length(pulses);

% construct parameter "key" to correspond to simulation numbers
% the order corresponds to loop order in multi_heat.sh (r-a-f-t)
T_key = zeros(rlength*plength*flength,3);
for ii = 1:rlength
   for jj = 1:flength
      for kk = 1:plength
         T_key((ii-1)*flength+(jj-1)*plength+(kk-1)+1,1) = radii(ii);
         T_key((ii-1)*flength+(jj-1)*plength+(kk-1)+1,2) = fluences(ii);
         T_key((ii-1)*flength+(jj-1)*plength+(kk-1)+1,3) = pulses(ii);
      end
   end
end

% Ep refers to induced particle electric field strength

%%%%%% Plotting of melt, eval, explosion %%%%%%%
Tmaxt=load('..//results//T_max.txt');  % fix this!
Tmaxe=Tmaxt(sim_min:sim_max,1);
Tmaxp=Tmaxt(sim_min:sim_max,2);
Tmaxw=Tmaxt(sim_min:sim_max,3);

epsilon=3*G./(g.*T_key(sim_min:sim_max,1));
beta=T_key(sim_min:sim_max,2).*T_key(sim_min:sim_max,4)./...
(Ceo*To*(4/3)*pi.*T_key(sim_min:sim_max,1).^3);
gamma=T_key(sim_min:sim_max,3)*g/Ceo;

neps=1;
leps=sim_tot/neps;

count=1;
dbeta=diff(beta);
for tt=1:leps-1
   if dbeta(tt)>0
      count=count+1;
   end
end
nbeta=count;
lbeta=leps/nbeta;
ngamma=lbeta;

vgamma=gamma(1:ngamma);
vbeta=zeros(nbeta,1);
for kk=1:nbeta
   vbeta(kk,1)=beta(kk*lbeta-1,1);
end

count=0;
Emat=zeros(ngamma,nbeta);
Emat_sc=zeros(ngamma,nbeta);
Tp=Emat;
Tw=Emat;
for oo=1:nbeta
   count=count+1;
   olow=(count-1)*ngamma+1;
   ohigh=(count)*ngamma;
   Emat(1:ngamma,oo,1)=Ep(olow:ohigh,1);
   Emat_sc(1:ngamma,oo,1)=Ep_sc(olow:ohigh,1);
   Tp(1:ngamma,oo,1)=Tmaxp(olow:ohigh,1);
   Tw(1:ngamma,oo,1)=Tmaxw(olow:ohigh,1);
end

Ecrit=2.7e10;%                    Critical efield Au monolayer ablation
Etot=Ecrit*ro/(latspace);%        Total particle ablation
Eone=ro*ec/(3*eo*4/3*pi*ro^3);%   Single electron emission
Tp_melt=1337;%                    Melt temp Au
Tw_bin=373;%                      Reg boil temp H2O (binodial)
Tw_spin=467;%                     Exp boil temp H2O (spinodial)

refine=5;
x1=(1:1/refine:length(vbeta));
lvbeta=interp1(vbeta,x1,'spline');
x2=(1:1/refine:length(vgamma));
lvgamma=interp1(vgamma,x2,'spline');
[xs,ys]=meshgrid(1:length(vbeta),1:length(vgamma));
[xb,yb]=meshgrid(x1,x2);
Ematb=exp(interp2(xs,ys,log(Emat),xb,yb,'cubic'));    
Twb=exp(interp2(xs,ys,log(Tw),xb,yb,'cubic')); 
Tpb=exp(interp2(xs,ys,log(Tp),xb,yb,'cubic')); 

figure(11)
contour(lvbeta,lvgamma,Ematb,[Ecrit Ecrit],'LineWidth',1,'LineColor','k');
xlabel('\beta','Fontsize',16);%,'FontWeight','BOLD');
ylabel('\gamma','Fontsize',16);%,'FontWeight','BOLD');
set(gca, 'xscale', 'log','yscale', 'log')
xlim([10^1 10^5]);
hold on

contour(vbeta,vgamma,Emat_sc,[Ecrit Ecrit]);
contour(lvbeta,lvgamma,Ematb,[Etot Etot],'LineWidth',1,'LineColor','k');
%contour(lvbeta,lvgamma,Twb,[Tw_bin Tw_bin],'LineWidth',1,'LineColor','k');
contour(lvbeta,lvgamma,Twb,[Tw_spin Tw_spin],'LineWidth',1,'LineColor','k');
%contour(vbeta,vgamma,Tw,[0.9*Tw_spin 0.9*Tw_spin]);
contour(lvbeta,lvgamma,Tpb,[Tp_melt Tp_melt],'LineWidth',1,'LineColor','k');
%legend('monolayer ablation','complete ablation',...
%         'water binodal','water spinodal','particle melting','Location','NorthWest');

%plot(1142,0.191,'ro');

figure(12)
lvf=lvbeta.*(Ceo*To*4/3*pi*(ro)^3/(Qabs*pi*(ro)^2));
lvt=lvgamma.*(Ceo/g);
contour(lvf,lvt,Ematb,[Ecrit Ecrit],'LineWidth',1,'LineColor','k');
xlabel('Fluence (J/m^2)','Fontsize',16);%,'FontWeight','BOLD');
ylabel('Pulse Duration (s)','Fontsize',16);%,'FontWeight','BOLD');
set(gca, 'xscale', 'log','yscale', 'log')
xlim([10^1 10^5]);
hold on

%contour(vbeta,vgamma,Emat_sc,[Ecrit Ecrit]);
contour(lvf,lvt,Ematb,[Eone Eone],'LineWidth',1,'LineColor','k');
contour(lvf,lvt,Ematb,[Etot Etot],'LineWidth',1,'LineColor','k');
contour(lvf,lvt,Twb,[Tw_bin Tw_bin],'LineWidth',1,'LineColor','k');
contour(lvf,lvt,Twb,[Tw_spin Tw_spin],'LineWidth',1,'LineColor','k');
%contour(vbeta,vgamma,Tw,[0.9*Tw_spin 0.9*Tw_spin]);
%contour(lvf,lvt,Tpb,[Tp_melt Tp_melt],'LineWidth',1,'LineColor','k');
%     legend('monolayer ablation','complete ablation',...
%         'water binodal','water spinodal','particle melting','Location','NorthWest');

%    plot(Fnp,taup,'bo','MarkerSize',7,'LineWidth',2);
%    plot(40,taup,'ro','MarkerSize',7,'LineWidth',2);
%herrorbar(Fnp,taup,2*Fnpe); 


